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We perform classical molecular dynamics to investigate the effects of mechanical strain on single¬ 
layer black phosphorus nanoresonators at different temperatures. We find that the resonant fre¬ 
quency is highly anisotropic in black phosphorus due to its intrinsic puckered configuration, and 
that the quality factor in the armchair direction is higher than in the zigzag direction at room 
temperature. The quality factors are also found to be intrinsically larger than graphene and M 0 S 2 
nanoresonators. The quality factors can be increased by more than a factor of two by applying ten¬ 
sile strain, with uniaxial strain in the armchair direction being most effective. However, there is an 
upper bound for the quality factor increase due to nonlinear effects at large strains, after which the 
quality factor decreases. The tension induced nonlinear effect is stronger along the zigzag direction, 
resulting in a smaller maximum strain for quality factor enhancement. 
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Black phosphorus (BP) is a new two-dimensional nano¬ 
material that is comprised of atomic layers of phosphorus 
stacked via van der Waals forces^. BP brings a number 
of unique properties unavailable in other two-dimensional 
crystals material. For example, BP has anisotropic prop¬ 
erties due to its puckered configuration^^”— 

While most existing experiments have been focused on 
potential electronic applications of BP—”—, a recent ex¬ 
periment showed that the resonant vibration response 
of BP resonators (BPR) can be achieved at a very high 
frequency^^ However, there have been no theoretical stud¬ 
ies on the intrinsic dissipation in BPRs to-date. In partic¬ 
ular, it is interesting and important to characterize the 
effects of mechanical strain on the quality (Q)-factors 
of BPRs given its anisotropic crystal structure, and fur¬ 
thermore considering that mechanical strain can act as 
an efficient tool to manipulate various physical proper¬ 
ties in the BP structureFor example, a large uni¬ 
axial strain in the direction normal to the SLBP plane 
can even induce a semiconductor-metal transition^i^”— 
We thus investigate the mechanical strain effect on the 
BPRs of armchair and zigzag directions, at different tem¬ 
peratures. 

In this work, we examine the effect of mechanical ten¬ 
sion on single-layer BPR (SLBPR) via classical molecular 
dynamical (MD) simulations. Both uniaxial and biaxial 
tension are found to increase the quality factor of the 
SLBPR, as the resonant frequency is enhanced by the 
applied tension. However, the Q-factor decreases beyond 


a critical strain value due to the introduction of nonlinear 
energy dissipation, which becomes dominant at large ten¬ 
sile strains. As a result, there is a critical strain at which 
the quality factor reaches the maximum value, which is 
about 4% and 8% at 50 K for mechanical tension along 
the zigzag and armchair directions, respectively. We find 
that the nonlinear dissipation is stronger if the BPR is 
stretched along the zigzag direction, which results in a 
smaller critical strain. 

Fig.[T]shows the structure of SLBP of dimension 50 x 50 
A that is used in our simulations. The atomic inter¬ 
actions are described by a recently-developed Stillinger- 
Weber potential^^^ The BPR simulations are performed 
in the following manner. First, the entire system is ther- 
malized to a constant temperature within the NPT (i.e., 
the particles number N, the pressure P and the tempera¬ 
ture T of the system are constant) ensemble by the Nose- 
Hoover— thermostat, which is run for 200 ps. Second, 
the SLBP is stretched by uniaxial or biaxial strain along 
the armchair or zigzag directions. The mechanical strain 
is applied at a strain rate of e = 0.0001 ps“^, which 
is a typical value in MD simulations. Third, the res¬ 
onant oscillation of the SLBP is actuated by adding a 
sine-shaped velocity distribution, uq sin(7rXj/L), to the 
system. In all simulations, we apply the velocity am¬ 
plitude uo = 2.0 A/ps, which is small enough to keep 
the resonant oscillation in the linear region. Fourth, the 
resonant oscillation of the SLBP is simulated within the 
NVE (i.e., the particles number N, the volume V and the 
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FIG. 1: (Color online) Configuration of SLBP with dimen¬ 
sions 50 X 50 A, from the top view in the top panel, and from 
the side view in the bottom panel. The total number of atoms 
is 660. 


temperature T of the system are constant) ensemble for 
90 ns, and the oscillation energy is recorded to extract 
the Q-factor. 

We first examine the intrinsic energy dissipation of the 
SLBPRs along the armchair and zigzag directions. The 
intrinsic energy dissipation is induced by thermal vibra¬ 
tions at finite temperatures. Fig. [2] shows the kinetic en¬ 
ergy time history in armchair SLBP at 4.2 K, 30 K and 
50 K. The oscillation amplitude of the kinetic energy de¬ 
cays gradually, which reflects the energy dissipation dur¬ 
ing the resonant oscillation of the SLBPR. As the temper¬ 
ature increases, the energy dissipation becomes stronger, 
indicating a lower Q-factor at higher temperature. 

The frequency and the Q-factor of the resonant oscilla¬ 
tion can be extracted from the kinetic energy time history 
by fitting to the function Ek{t) = + cos(27r2/t)(l — 

• The first term represents the averaged kinetic 
energy after the resonant oscillation has completely de¬ 
cayed. The constant E^ is the total kinetic energy at 
t = 0, i.e. at the moment when the resonant oscilla¬ 
tion is actuated. The frequency of the resonant oscilla¬ 
tion is /, so the frequency of the kinetic energy is 2/. 
The kinetic energy time history is usually a very long 
data set, so it is almost impossible to fit it directly to 
the above function. The fitting procedure is thus done 
in the following two steps as shown in Fig. [3l First, 
Fig. [3] (a) shows that the energy time history is fitted 
to the function Ek{t) = E/^ + E^ cos(27r2/t) in a very 
small time region t G [0, 50] ps, where the approximation 
(1 — ^ 1 has been done for the Q-factor term as the 

energy dissipation is negligible in the small time range. 
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FIG. 2: (Color online) The kinetic energy per atom for arm¬ 
chair SLBPR at 4.2 K, 30 K and 50 K from top to bottom. 
The Q-factors are 1621900, 110000 and 63250 respectively. 


The parameters E/^, E^, and / are obtained accurately 
from this step. Second, Fig. [3](b) shows that the oscilla¬ 
tion amplitude of the kinetic energy can be fitted to the 
function E™^(t) = E^(l — ^)'^Mn the whole simulation 
range t G [0, 90] ns, which determines the Q-factor. Fol¬ 
lowing these fitting procedures, the Q-factor is 63250 for 
the armchair SLBPR at 50 K. 

Fig. 0] shows the temperature dependence for the Q- 
factor of the SLBPR along the armchair and zigzag direc¬ 
tions. At most temperatures, the Q-factor is larger in the 
armchair direction. It means that the energy dissipation 
is weaker for armchair SLBPR, considering that the fre¬ 
quency in the armchair SLBPR is only half of that in the 
zigzag SLBPR^2^ The temperature dependence for the 
Q-factor can be fitted to the function Q = 1.9 x 
and Q = 3.0 x for armchair and zigzag SLBPR, 

respectively. 

These Q-factors are higher than the Q-factors in 
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FIG. 3: (Color online) Two-step fitting procedure to extract 
the frequency and Q-factor from the kinetic energy time his¬ 
tory for armchair SLBPR at 50 K. (a) The kinetic energy is 
htted to the function Ek{t) — Ek + cos(27r2/t) in a small 

time range, giving the frequency / = 0.090874 THz. (b) The 
kinetic energy is htted to the function = E^{1—^)^^ 

in the whole time range, yielding the Q-factor value of 63250. 


graphene nanoresonators {Q = 7.8 x iQ4y-i-2) ^25^26 xhis 
is likely because there is also a large energy bandgap in 
the phonon dispersion of SLBP^ which helps to pre¬ 
serve the resonant oscillation of the SLBPR^^^ In con¬ 
trast, there is no such energy bandgap in the phonon 
dispersion of graphene, so the SLBPR has higher Q- 
factor than graphene nanoresonators. The Q-factors of 
SLBPR’s are also higher than those of M 0 S 2 nanores¬ 
onators (Q = 5.7 X Both SLBP and M 0 S 2 

have energy bandgaps in their phonon dispersions. This 
is important as our simulation results imply that non¬ 
linear phonon-phonon scattering is weaker in SLBP, i.e., 
the resonant energy dissipation is weaker in SLBP than 
M 0 S 2 . 

We now report the effects of mechanical strain on 
both armchair and zigzag SLBPRs at 50 K. We consider 
four cases, i.e., (I) the effect of uniaxial strain on arm¬ 
chair SLBPR, (II) the effect of uniaxial strain on zigzag 
SLBPR, (III) the effect of biaxial strain on armchair 
SLBPR, and (IV) the effect of biaxial strain on zigzag 



FIG. 4: (Golor online) Temperature dependence for the Q- 
factors of armchair and zigzag SLBPRs. 



FIG. 5: (Golor online) Strain dependence for the Q-factor 
of SLBPR in four cases at 50K. The Q-factor depends on 
strain as the function Q/Qo = —ae^ + 5e + 1.0, which gives a 
maximum Q-factor value at a critical strain. 


SLBPR. Fig. [ 5 ] shows the strain dependence for the Q- 
factor (with reference to the value Qo without strain) of 
SLBPR under uniaxial or biaxial mechanical tension. In 
case I, the mechanical strain is applied purely in the arm¬ 
chair direction, while the SLBP is stretched in the zigzag 
direction in the other three cases. 

For all of the four cases, the Q-factor first increases and 
then decreases after a critical strain value. The Q-factor 
depends on the strain as the function Q/Qo = — ae^ + 6e-f 
1.0, where the fitting parameters (a, 6) are (0.029, 0.42), 
(0.055, 0.40), (0.043, 0.36), and (0.070, 0.63) for the four 
studied cases, respectively. The linear term he represents 
the enhancement effect on the Q-factor by the mechanical 
tension, as the frequency of the resonator is increased by 
the tension in the small strain range. The quardratic 
term —ae^ is because the Q-factor will be reduced by 
the nonlinear effect resulting from the mechanical tension 
in the large strain range. The interplay between these 
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FIG. 6: (Color online) Stress-strain relation for SLBP under 
mechanical tension. The stress (cr) is fitted to a function of 
strain (e) as cr = Ee + ^ with E as the Young’s modulus 

and D as the TOEC. The nonlinear effect is estimated by the 

ratio 7 = 



EIG. 7: (Color online) Strain dependence for the Q-factor of 
SLBPR in four cases at 170K. 

two competing effects results in a maximum value for 
the Q-factor at a critical strain Cc- The critical strain 
value is about 8 % for case I, in which the mechanical 
tension is applied only in the armchair direction. For all 
other three cases, the critical strain is around 4%, where 
the mechanical tension has a component in the zigzag 
direction. 

The differences in the above critical strains can be 
understood from the different strain induced nonlinear 


properties in the SLBP. Fig. [ 6 ] shows the stress-strain 
curve for the SLBP stretched in the above four cases. 
The stress-strain curve is fitted to the function a = 
Es + \De‘^ ^ with E as the Young’s modulus and D as 
the third-order elastic constant (TOEC)^. The nonlin- 

-D 

ear to linear ratio of 7 = gives an overall estimation 
of the strain induced nonlinear effect on the SLBP. The 
parameter 7 is found to be -1.66 for case I, -3.64 for case 
II and -3.65 for the other two cases. This means that 
the nonlinear effect is the weakest in case I, where the 
SLBP is stretched purely in the armchair direction. As a 
result, the parameter a has the smallest value for case I, 
leading to the largest critical strain. This phenomenon (a 
maximum Q factor due to the strain effect) has also been 
obtained in nanowire resonators. For example, Kim and 
Park found that the maximum Q factor occurs around 
1.5% tensile strain in the metal nanowire resonators^^ 

Fig.[7]shows the strain effect on the Q-factor at 170 K 
for all four cases. The critical strain is also observed at 
this higher temperature, and the critical strain value for 
SLBR at 170K is about 5% for case I and around 2-3% 
for other three cases. However, the critical strain value 
is smaller as compared with the critical strain at 50 K in 
Fig. [5l This is because the nonlinear effect is stronger at 
higher temperature due to the thermally-induced random 
vibrations. The combination of the two nonlinear effects 
(induced by temperature and strain) leads to a smaller 
critical strain at higher temperature. 

In conclusion, we have performed classical molecular 
dynamics simulations to study the effects of mechanical 
tension effects on the SLBPR at different temperatures. 
We find that intrinsically, or neglecting strain, the Q- 
factors for armchair SLBPR are generally higher than 
for zigzag SLBPR, and are also larger than those found 
previously in graphene and M 0 S 2 nanoresonators. When 
the effects of mechanical strain are considered, our key 
finding is that there is a maximum point in the strain de¬ 
pendence of the Q-factor due to the competition between 
the enhancement at small strains and the reduction due 
to nonlinear effects at large strains. 
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